TCGA BRCA

library(TCGAbiolinks)

# query <- GDCquery(project = "TCGA-BRCA",
#                   data.category = "Transcriptome Profiling",
#                   data.type = "Gene Expression Quantification")
# 
query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  sample.type = c("Primary Tumor","Solid Tissue Normal"))

df = query$results[[1]]
df$sample_type %>% table() %>% data.frame()

# Download and prepare data
GDCdownload(query)
data <- GDCprepare(query)

library(DESeq2)

# 가능한 모든 assays 이름을 확인
assayNames(data)

# unstranded: Strand 비특이적 발현 데이터
# stranded_first: Strand 특이적 (첫 번째 strand) 발현 데이터
# stranded_second: Strand 특이적 (두 번째 strand) 발현 데이터
# tpm_unstrand: TPM (Transcripts Per Million) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_unstrand: FPKM (Fragments Per Kilobase of transcript per Million mapped reads) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_uq_unstrand: FPKM-UQ (Upper Quartile normalized FPKM) 방식으로 정규화된 strand 비특이적 발현 데이터

# 'unstranded' 데이터를 사용하여 발현 수치 추출
counts <- assay(data, "unstranded")

# ENSG to Symbols
library(org.Hs.eg.db)

# 유전자 ID에서 버전 제거
ensembl_ids <- gsub("\\..*", "", rownames(counts))

# 중복된 ID 처리
unique_ids <- unique(ensembl_ids)

# 중복 제거된 ID에 대해 유전자 심볼 매핑
gene_symbols <- mapIds(org.Hs.eg.db,
                       keys = unique_ids,
                       column = "SYMBOL",
                       keytype = "ENSEMBL",
                       multiVals = "first")

# 원래 데이터에 매핑된 유전자 심볼 할당
# 여기서 중복 제거된 목록을 이용해 각각의 원본 ID에 매핑된 심볼 할당
symbol_names <- gene_symbols[ensembl_ids]
rownames(counts) <- symbol_names

# Metadata 
filtered_colData <- colData(data)[colnames(counts),]

TCGA_BRCA_countsMeta = list(counts = counts,
                            meta = filtered_colData)
TCGA_BRCA_countsMeta %>% saveRDS(paste0(dir,"TCGA_BRCA_countsMeta.rds"))

# Generate info table
info <- data.frame(matrix(nrow = ncol(count.mtx), ncol = 2))
colnames(info) <- c('sample', 'cond')
info$sample <- colnames(count.mtx)
info$cond <- dds$sample_type
info$cond <- factor(info$cond, 
                    levels = c("Solid Tissue Normal","Primary Tumor")) # CTL going first
# levels(info$cond)

# DESeq
dds <- DESeqDataSetFromMatrix(count.mtx, info, ~ cond)
dds <- DESeq(dds) 
# dds %>% saveRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)

dds = readRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)

# Add DEG information 
fc = 1.5
pval = 0.05

# res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & pvalue < pval, 'UP',
#                                ifelse(log2FoldChange <= -log2(fc) & pvalue < pval, 'DN','no_sig')))

res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & padj < pval, 'UP',
                               ifelse(log2FoldChange <= -log2(fc) & padj < pval, 'DN','no_sig')))
res = na.omit(res)

res %>% DT::datatable()

res$DE = factor(res$DE, levels = c('UP','DN','no_sig'))
res %>% 
  ggplot(aes(log2FoldChange, -log10(padj), color=DE)) + 
  geom_point(size=1, alpha=0.5) + 
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_classic() +
  geom_vline(xintercept = c(-log2(fc),log2(fc)), color='grey') +
  geom_hline(yintercept = -log10(0.05),color='grey') +
  guides(colour = guide_legend(override.aes = list(size=5))) +
  ggtitle(paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )) +
  ggeasy::easy_center_title() ## to center title

t= paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )
up = nrow(res[res$DE == "UP", ])
dn = nrow(res[res$DE == "DN", ])
res %>% ggplot(aes(log2FoldChange, -log10(padj), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(fc), log2(fc)), size=0.1, color="grey") +
  geom_hline(yintercept = -log10(0.05), size=0.1, color="grey") +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

res %>% filter(DE != "no_sig") %>% 
  ggplot(aes(DE, fill=DE)) + geom_bar(color="black", size=0.2) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size= 4, color=c("salmon","royalblue")) +
  scale_fill_manual(values = c("salmon", "royalblue"), guide=F) +
  theme_bw()

library(clusterProfiler)
hallmark <- msigdbr::msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol)

perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$log2FoldChange
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  result <- x@result %>% arrange(desc(NES))
  result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
  return(result)
}

# Application 
gsea.res = perform_GSEA(res = res, ref = hallmark) 

filtered_gsea = gsea.res %>% mutate(sig= ifelse(pvalue <= 0.05,"p value <= 0.05", "p value > 0.05"))

# Modified GSEA NES plot 
gsea_nes_plot =function(gsea.res, title, fontsize.x = 5, fontsize.y = 6){
  gsea.res %>% ggplot(aes(reorder(ID, NES), NES)) +
    geom_col(aes(fill=sig), color="grey1", size=0.2) +
    coord_flip() +
    labs(x="Pathway", y="Normalized Enrichment Score",
         title= "GSEA") + 
    theme_classic() +
    # scale_fill_gradient(low = '#FF0000', high = '#E5E7E9') +
    scale_fill_manual(values = c("#FF0000","grey88")) +
    theme(axis.text.x= element_text(size=fontsize.x, face = 'bold'),
          axis.text.y= element_text(size=fontsize.y, face = 'bold'), 
          axis.title =element_text(size=10)) +ggtitle(title)
}

t= paste0("DEGs from ",levels(dds$cond)[2], " / ", levels(dds$cond)[1] )
gsea_nes_plot(filtered_gsea, title = t, fontsize.x = 8, fontsize.y = 8)

gsea.res %>% DT::datatable()